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ABSTRACT 

We present the results of a systematic numerical study of the onset of mass transfer in double 
degenerate binary systems and its impact on the subsequent evolution. All investigated systems be- 
long to the regime of direct impact, unstable mass transfer. In all of the investigated cases, even 
those considered unstable by conventional stability analysis, we find a long-lived mass transfer phase 
continuing for as many as several dozen orbital periods. This settles a recent debate sparked by a 
discrepancy between earlier SPH calculations that showed disruptions after a few orbital periods and 
newer grid-based studies in which mass transfer continued for tens of orbits. The number of orbits 
a binary survives sensitively depends on the exact initial conditions. We find that the approximate 
initial conditions that have been used in most previous SPH calculations have a serious impact on all 
stages of the evolution from the onset of mass transfer up to the final structure of the remnant. We 
compare "approximate" initial conditions where spherical stars are placed at an initial separation ob- 
tained from an estimate of the Roche lobe size with "accurate" initial conditions that were constructed 
by carefully driving the binary system to equilibrium by a relaxation scheme. Simulations that use the 
approximate initial conditions underestimate the initial separation when mass transfer sets in, which 
yields a binary that only survives for only few orbits and thus a rapidly fading gravitational wave 
signal. Conversely, the accurate initial conditions produce a binary system in which the mass transfer 
phase is extended by almost two orders of magnitude in time, resulting in a gravitational wave signal 
with amplitude and frequency that remain essentially constant up until merger. The mass transfer 
also shows a unique oscillatory signal that is most pronounced for binary systems that are marginally 
unstable. As we show that these binaries can survive at small separation for hundreds of orbital peri- 
ods, their associated gravitational wave signal should be included when calculating the gravitational 
wave foreground (although expected to below EISA's sensitivity at these high frequencies). We also 
show that the inclusion of the entropy increase associated with shock-heating of the accreted material 
reduces the number of orbits a binary survives given the same initial conditions, although the effect 
is not as pronounced when using the appropriate initial conditions. The use of accurate initial condi- 
tions and a correct treatment of shock heating allows for a reliable time evolution of the temperature, 
density, and angular momentum, which are important when considering thermonuclear events that 
may occur during the mass transfer phase and/or after merger. Our treatment allows us to accurately 
identify when surface detonations may occur in the lead-up to the merger, as well as the properties 
of final merger products. 

Subject headings: supernova: general - WDs - nuclear reactions, nucleosynthesis, abundances - hy- 
drodynamics 



1. INTRODUCTION 



White dwarfs (WDs) are the end point of the evolution of the vast majority of the sta rs in the universe. Our Galaxy 
harbours of order 10^° WDs (Napiwot zkill2QQ9|) and about 2.5 x 10^ close binary WDs (Nelemans et al.||2QQlb ), often 
referred to as double dee;enerates (DD). For about half of them a Hubble time is Ions; enoue;h so that gravitational 



{DDJ 

wave emission can drive th em to a phase of mass transfer. Such interacting DDs are suspected to produc e observable 
objects such a s sdb stars (|Webbink||1984| llben fc Tutukov||1986| |Saio fc Jeffery"2QQQ> W n et al.|[2Q^, R Corona 
Borealis stars ( |Webbink|[l984[ [Clayton et al. 2007 K AM CVn stars (Paczyhski 1967; Nelem anTet al.||2U(Jlal [Solheim 
2010 ) and, maybe most spectacularly, type la supernovae ( Iben fc TutukoV|1984^|Webbink|1984| ). While tor many years 
there was a strong in clination towards the single dege nerate scenario, the recent disc overy of super-Chandrasekhar" 
events (Howe ll et al.||2Q06| |Hicken e t al. 2007| |Yuan et al.|2007 Silverman et al.|2011 ) has re-strengthened the interest 
in DD mergers as possible type la progenitors! I'he lead- up to a merger includes a phase of mass transfer that is 
crucial for understanding the evolution and fates of double WD systems. This phase will affect the orbital evolution, 
the angular momentum and thermodynamical state of the final merger product, the associated gravitational wave 
signal, the environment into which any potential ejecta produced by an thermonuclear explosion will expand into, and 
whether or not a particular binary will merge at all. The numerical study of the onset of this mass transfer phase is 
the main topic of this investigation. 
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rosswog@j acobs-uni versity. de 

^ TASC, Department of Astronomy and Astrophysics, University of Cahfornia, Santa Cruz, CA 95064; jfg@ucohck.org, en- 
rico@ucohck.org 
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The merger process of two WDs has been modeled by a number of groups ^Benz et al. ||199Q[ |Rasio fc Shapiro||1995l 
Segretain et al. 1997 Guerrero et al.||2004| |Yoon et al.||2007| |Pakmor et al. 2010), all of which use the smoothec. 



particle hydrodynamics (SFH) method. In all of these s imulati ons the mass transfer phase lasted for very few orbital 



time scales. Recently, grid-based simulations ([Motl et al. 2002; D'Souza et al.| [2QQ6{ Motl et aL]|2007[ ) were performed 
in which the authors carefully tried to reduce the angular momentum non-conservation due to advection errors that 
often plague grid-based codes ( New fc Tohline||1997 ). The latter results differed from the previous SPH calculations 
in that they showed long-lived mass transfer for many orbital periods and this has sparked some discussion about the 
stability of mass transfer in double degenerate sy stems during these late stages. Based on our experience in modeling 
the dynamics of neutron star black hole systems ( Rosswog et al.|[2QQ4 Rosswog||2005 ), we had suspected that double 
WD binaries may be very sensitive to the exact initial conditions. As we will demonstrate below, the initial conditions 
used have a decisive impact and are one of the major reasons for the disagreement between various mass transfer 
calculations. 
This is the companion paper of 



unstable, direct impact systems. 



Guillochon et al.| (I2010D 
he work showed t hi 



in which we focused on the fate of accreted helium in 
lat, contrary to prior belief, the mass transfer phase between 
two WDs may be eventful, and can be accompanied by Kelvin-Helmholtz-triggered thermonuclear surface explosions. 
These explosions may lead to the ignition of the accreting CO- WD in a double detonation scenario, which could be 
the long sought-after mechanism to initiate a type la supernova from WD binary merger. For the described work we 
had chosen a hybrid approach by combining result s using the SPH method for the orbital evolution and the adaptive 
mesh refinement code FLASH ( Fryxell et al.|[2QQQ ) to focus on the evolution of the low-density helium material. We 
consider this a "best-of-both-worlds" approach, because it combines the major strengths of both Lagrangian and 
Eulerian methods. In this paper we present the details of the mass transfer/orbital evolution calculations that have 
been obtained with SPH, as well as a detailed comparison with the previously published FLASH calculations. 

The main focus of this paper is the onset of the mass transfer and its impact on the orbital evolution; the structure 
of the resulting remnants will be discussed elsewhere. In Section [2] we briefiy summarize conditions for the stability of 
mass transfer, in Section |3] we focus on the numerical aspects of the mass transfer modeling. We describe in detail how 
we construct our initial conditions and demonstrate their impact on the outcome of essentially every aspect of double 
degenerate merger simulations. We also briefiy summarize the essential features of the SPH code that we are using. In 
a set of test calculations we explore to which extent the numerical resolution, the form of artificial viscosity, and the 
exact numerical initial conditions affect the results. In Section [4] we highlight results from two systems representative 
of a set of binary evolution calculations. In Sections [5] and [6] we discuss the implications of the extended mass transfer 
phase on the gravitational wave signal and the triggering of surface detonations prior to merger, respectively. We 
conclude the paper in Section [7| with a summary. 



2. STABILITY OF MASS TRANSFER 

In what follows we briefiy collect the standard arguments on the stability of mass transfer (see e.g. 'Fran k et al 
2002 ) . Once stars are close enough that mass transfer via Roche lobe overfiow can set in, the further evolution of the 



binary system depends on how both the orbit/Roche lobe size and the star react to the re- arrangement of mass and 
angular momentum. Due to their particular mass radius relation, WDs react to mass loss by expansion which tends 
to enhance the mass transfer rate. If mass is transferred to the heavier star, the orbit has to widen to ensure the 
conservation of the centre of mass and, therefore, mass transfer is tendentially quenched. The interplay of both effects 
determines the detailed orbital evolution. The total angular momentum of a point mass binary is given by 



J = M1M2 



Ga 



1/2 



(1) 



where the Mi are the masses, G is the gravitational constant, Mtot is the total mass and a the separation of the binary 
system. In the remainder of this work we will always denote the donor by label 2, the accretor by 1. On logarithmic 
differentiation one finds 



a 

2a 



J 
J 



M2 



(2) 



where q = M2/M1. As expected, angular momentum loss tends (J < 0) to shrink the orbit, while mass transfer 
(M2 < 0) tends to increase the orbital separation. The change of the donor Roche lobe size c an be related to the 
change in orbital separation by logarithmically differentiating a simple Roche lobe size estimate ( Paczyhs"ki||1971 ) 



Rl,2 = 0.462 a 



1/3 



(3) 



which yields 



L,2 



Rl,2 



IM2 
3M2' 



(4) 
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Using this to eliminate a from Equation ([2| then yields 

i^L,2 _ 2 j 2M2 ( 5 



i^L,2 J M2 y 6 



(5) 



and therefore for a mass ratio g > 5/6 a non-zero mass loss will trigger an expansion of the donor WD and at the same 
time shrink the Roche lobe size i^L,2- This will be further accelerated by any additional angular momentum loss, so 
with the used assumptions the mass transfer is clearly unstable in such cases. 
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Fig. 1. — Analytical estimates of mass transfer stability in double WD binaries. The region below the lower solid line indicates the 
regime where a disk forms. Above the line, the accretion stream directly impacts on the accretor. The dot-dashed line separates systems 
that undergo sub- vs. super-Eddington accretion, with super-Eddington accretion possibly leading to an unstable binary. The hatched 
region of the plot indicates systems for which the accretion is both sub-Eddington and through a disk, which are almost certainly stable. 
For mass ratios q larger than 2/3, the mass transfer is guaranteed to be unstable. The area between the regions of guaranteed stability 
and instability corresponds to the systems that can be either stable or unstable , depending on the spin-orbit coupling. The filled squares 
indicate the masses of the systems simulated in this work, seefl] Adapted from Marsh et aL](|2004|). 



Of course, the so far discussed case is highly idealized. A point mass binary was assumed where the mass lost by the 
donor was completely incorporated into the accretor, while in reality mass can be lost from the system. Depending on 
whether or not the circularization radius exceeds the accretor radius, the transferred mass can either form a disk or 
it can impact directly onto the accretor surface. If a disk is formed, it would be tidally perturbed by the companion 
star and can thus feed back some fraction of the angular momentum into the orbit. In the direct impact case most 
of the angular momentum is invested into spinning up the accretor, which removes orbital angular momentum from 
and is therefore expected to destabilize the orbit. It should be noted, however, that a deformed donor can also feed 
back some angular momentum into the orbit. The orbit and the diffe rent binary components can exchange angular 
momentum both by advection of mass and by tides. Marsh et al. (2004) have generalized the simple treatment leading 
to Equation ([2| for direct impact and tidal spin evolution 
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2a 
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Jorh 



J1 + J2 



orb 



tid 



M2 
M2 



1-q - ^{l+q)r. 



(6) 



where Vc is the circularization radius in units of the orbital separation. The second term accounts for the exchange 
of angular momentum due to tides, the third term accounts for the spin of the accretor. If tides and accretor spin 
are ignored. Equation (|6| reduces to In Figure [l] we show the stability criteria for binaries o f different donor and 
accretor masses, superimposed with the initial masses of our simulations (Table [T]). [Marsh et al. ( 2004 J find the upper 
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left of the diagram to be guaranteed unstable, the lower right guaranteed stable and an uncertain region along the 
diagonal of the diagram. In this work we mainly investigate the "uncertain" region in the diagram. 

3. NUMERICAL MODELING OF MASS TRANSFER 

The main tool of this investigation is a 3D sm oothed particle hydrodynamics (SPH) code. For recent reviews of the 
method see e.g. Monaghan (2005) and Rosswog (2009). The orbital dynamics of a binary system is very sensitive to 
the (re-) distribution of angular momentum and therefore a numerical simulation has to ensure very accurate angular 
momentum conservation. Its p urely La^ranffla n nature and the built-in conservation of mass, energy, linear and angular 
momentum (e.g. Section 2.4 in Rosswog|2009 ) make SPH a very powerful tool for this type of study. The conservation 
is exact up to effects that result from finite accuracy in gravitational forces due to the use of a tree, see below, and 
in the numerical integration of the ODE set. Both issues, however, are fully under control and can be improved to 
a desired level of accuracy by choosing the corresponding parameters for the tree and the ODE integration in an 
appropriate way, at the expense of a larger computational effort, though. For our parameters, even in the 84-orbit 
PI, corresponding to as many as 16,930 dynamical time scales of the 0.8 Mq star, energy and angular momentum 
are conserved to better than 1%. It has to be stated, however, that the advantage of exact mass conservation via 
constant SPH particle masses turns into a disadvantage if we ar e interes ted in the resolution of low-density portions 
of the flow. That is why we turned to the AMR code FLASH in [Guillochon et al.| ( |2010p to investigate the settling of 
the transferred material onto the accretor. 



Details of our SPH code can be found in Rosswog et al. (2008). It uses an artificial viscosity scheme ( Morris 
[Monaghan* 1997|) that reduces the dissipation terms to a very low level away from discontinuities, together witn a switch 
(fiSalsara 1995) to suppress the spurio us viscous forces in pur e shear flows. The system of fluid equations is closed by 
the HELMHOLTZ equation of state ( Timmes fc Swesty|2000 ). It accepts an externally calculated nuclear composi tion 
and allows a convenient coupling to a nuclear reaction network. We use a minimal nuclear reaction network (Hix| 
et al.|p!998 ) to determine the evolution of the nuclear composition and to include the energetic feedback onto the gas 



from the nuclear reactions. A set of only seven abundance groups greatly reduces the computational burden, bu t still 
reproduces the energy generation of all burning stages from He burning to NSE accurately. We use a binary tree ( Benz 
et al.||199Q) to search for the neighbor particles and to calculate the gravitational forces. 



Run 
No. 


Masses 


Comp. 


[10^] 


0-0,9 


[s] 


^max,8 


Pmax,7 


[%I^tot] 


[lO-^Mo] 


iVorb 


a'dis,9 


Mtrans 

[Mq] 


Comment 


Production runs 


PI 


0.2 + 0.8 


He-CO 


2 


5.77 


239 


4.89 


0.95 


5.7 


8.6 


84 


6.15 


0.08 


See|3.2.5||4.l| 


P2 


0.3 + 1.1 


He-CO 


2 


4.75 


151 


12.23 


6.15 


7.11 


17.6 


59 


5.01 


0.13 




P3 


0.5 + 1.2 


He-ONeMg 


2 


3.35 


81 


18.85 


15.44 


4.86 


17.4 


31 


3.27 


0.09 




P4 


0.3 + 0.6 


He-CO 


2 


4.04 


148 


4.92 


0.32 


2.5 


3.7 


45 


4.07 


0.07 




P5 


0.6 + 0.9 


CO-CO 


2 


2.68 


62 


10.6 


1.75 


0.8 


1.85 


29 


2.58 


0.10 


See|3.2.5!|4.2|^ 


G3 


Guillochon et al?| 


P6 


0.2 + 0.3 


He-He 


2 


4.39 


224 


2.37 


0.05 


1.03 


0.05 


22 


4.29 


0.05 




P7 


0.3 + 0.4 


He-He 


2 


3.60 


141 


2.72 


0.11 


0.6 


0.65 


14 


3.41 


0.06 




P8 


0.9 + 1.2 


CO-CO 


2 


1.91 


31 


71.15 


16.1 


2.3 


34.1 


29 


1.82 


0.07 




Test runs 


Gl 


0.45 + 0.9 


He-CO 


1 


3.36 


91 


16.11 


1.59 


2.2 


4.6 


30 


3.35 


0.12 


Guillochon et al. 




G2 


0.45 + 0.67 


He-CO 


1 


3.12 


90 


12.3 


0.48 


0.55 


0.9 


21 


3.01 


0.10 


jGuillochon et al. 




Tl 


0.6 + 0.9 


CO-CO 


2 


2.31 


49 


16.4 


1.83 


0.48 


2.01 


1 


2.02 


0.15 


app. IC 


T2 


0.2 + 0.8 


He-CO 


2 


4.95 


190 


5.24 


0.96 


1.91 


1.8 


2 


4.43 


0.08 


app. IC 


T3 


0.3 + 0.6 


He-CO 


2 


3.89 


139 


5.59 


0.33 


0.4 


0.5 


6 


3.64 


0.08 


ao = 1.13a^ss 


T4 


0.3 + 0.6 


He-CO 


0.4 


4.00 


145 


4.99 


0.35 


2.1 


2.6 


20 


4.0 


0.07 




T5 


0.3 + 0.6 


He-CO 


0.4 


4.00 


145 


6.10 


0.36 


1.91 


2.5 


23 


3.99 


0.08 


AV2 


T6 


0.3 + 0.6 


He-CO 


0.4 


4.00 


145 


4.70 


0.33 


2.09 


2.6 


18 


4.00 


0.08 




T7 


0.3 + 0.6 


He-CO 


0.4 


4.00 


145 


3.99 


0.31 


2.13 


2.5 


19 


4.12 


0.07 





TABLE 1 

Summary of the runs performed for this paper. ao,9 is the initial separation in units of 10^ cm, rmax,8 is the peak 

TEMPERATURE IN UNITS OF 10^ K, Pmax,7 THE PEAK DENSITY IN UNITS OF lO'^ G CM~^ AND L^nb IS THE FRACTION OF ANGULAR 
MOMENTUM CONTAINED IN MATERIAL THAT IS UNBOUND AT THE END OF THE SIMULATION. adis,9 IS THE SEPARATION (iN 10^ CM) AT WHICH 
TIDAL DISRUPTION OF THE DONOR OCCURS AND Mtrans IS THE AMOUNT OF MATERIAL (iN Mq) LOST BY THE DONOR PRIOR TO DISRUPTION. 



3.1. Initial conditions 

In what follows we will summarize how initial conditions are constructed for double WD binary systems, both in 
this and in previous work. 

3.1.1. Previous work 

In their pioneering work, 'Benz et al. ( 1990") followed the dynamical evolution of a 0.9 + 1.2 Mq WD binary system. 



They constructed cold, isolated WDs m hydrostatic equilibrium and placed them without spin-rotation at an initial 
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Fig. 2. — Equilibrium configuration (0.9 + 1.2 Mq) at the onset of mass transfer. The figure shows a projection of the SPH particles onto 
the orbital plane. Each particle's location within this plane is shown as a translucent grey circle. 

separation such that the seconda ry's Roche radius was smaher than its unperturbed initial radius by 8%. One of the 
major conclusions of Benz et alJ was that within about two orbital periods the secondary is completely destroyed and 
transformed into a thick accretion disk around the accretor. As we will show below, this short time to merger is an 
artifact of the initial separation being too small. We have followed the recipe of |Benz et al.| for the initial conditions 
and find exactly the same result as presented in their paper. We have further compared the orbital evolution of a 
50,000 particle bina ry sys tem, once according to the original Benz et al. description and according to our scheme, as 



described in Section 3.1.2 Again, with their initial conditions the merger set in within about one orbital period. With 



our in itial conditions and e veryth ing else be i ng th e same, the binary survived for 14 orbits of continued mass transfer. 



In Segretain et al. 
fills its Roche lobe. . 



( |1997[) and|Yoon et al.| (2007) the stars were placed at a mutual separation so that the secondary 
n|!Segretain et al. | the wDs had n o initial spin, while in Yoon et al. the WDs begin tidally locked. 



Guerrero et al. (2004) and Loren-Aguilar et al. (2009) initially placed the stars in a detached configuration and then 
slowly decreased the separation to the point where mass transfer began. However, their relaxation procedure did not 
allow the stars to reach complete hydrostatic equilibrium. If the stars are brought together too quickly, this can cause 
undesired (and unrealistic) oscillations in the WDs and may again lead to a n abbreviat ed mass transfer phase. For 
example, the mass transfer phase in their 0.4 + 1.2 Mq system (Figure 10 in Guerrero et al. (2004)) only lasts for a 
single orbit . This is also the case for the simulation of a 0.6 + 0.8 Mq systeni(see~Figures 1 andTY in |Loren-Aguilar| 

etal][2QQ9] ). 

Our in itial conditions are constructed in a way that is similar to [Rasio fc Shapiro ( 1995[ ) (Section [3.1.2 ). Rasio 



Shapiro] constructed a synchronized binary system with mass ratio of g = 0.5, and also followed the dynamical evolution 



once a few particles had crossed the Li point. They found a merger after five orbital periods and also questioned the 
earlier use of non-equilibrium initial conditions. We have tried to reproduce their calculations as closely as possible 
by adopting the same particle number, equation of state, and the same prescription f or evolving the entrop ic function 
K{s) (Equation (12)), but still we only found a merger after fifteen orbital periods. Fryer & Diehl (2008) place the 
two stars close eno ugh such that the do nor is already overfilling its Roche radius, which leads to a premature merger. 
The recent work of Pakmor et al. (2010) focuses primarily on the coalescence of the binary, but their initial conditions 
are not appropriate lor studies of long-term stability. 

3.1.2. Our approach 



We primarily focus on synchronized binary syste ms, but we also 



)erfor m some tests with initially non-spinning WDs. 
2004|). We begin our initialization procedure with 
K) equilibrium models. These stars are then mapped 



Our procedure is similar to the one described in Rosswog et al. 
stars that are set up according to one-dimensional, cold (T = 10^ 
into three dimensions and are then driven to an accurate numerical equilibrium by applying a velocity-dependent 
acceleration /damp while keeping the temperature constant. The stars are then placed at an initial separation that is 
large enough to avoid any immediate mass transfer. The full SPH code is then used to relax the initial conditions to a 
tidally deformed binary system at the brink of numerically resolvable mass transfer. This phase of the initialization is 



6 



Dan, Rosswog, Guillochon and Ramirez-Ruiz 




x[W cm] 

Fig. 3. — Projection of the SPH particles on the (x, <l>pm/^) plane. The dotted red line is the point mass Roche potential, ^pm(x,y = 
O^z = 0). The solid black line is the value of the potential ^ along the x— axis, the SPH particle values are shown as filled black circles. 
Left: Wide view showing both stars. Right Zoom into the neighborhood of the Lagrange point Li. The fluid potential ^ lies below the 
poiivt massivalue ^urn and .therefore the point niass. approxirnatipn .underestimates the -distance where) maffltr^nsfer setSiin. .i i i 

perlormed m a corotatmg frame, i.^. centrimgai and Corions lorces must De accounted lof^l We calculate the orbital 




frequency cj in a way that the total forces on the individual stellar centers of mass vanish exactly (Rosswog et aH2Q04 ) 
and thus account for finite size deviations from the Kepler frequency. The acceleration applied to a particle t then 
reads 

fi = /grav,i + /hyd,i - X [iV X n) - 2u; X Vi . (7) 



Here, /grav,i and /hyd,i are the gravitational and hydrodynamic acceleration, Vi and Vi position and velocity in the 
corotating frame, fli^2 is either the primary or the secondary star's angular velocity, and a damping time scale. In 
this frame, fti^2 = for synchronized systems and fti^2 = for non-spinning stars. During this relaxation process, 
we adiabatically decrease the orbital separation a in a way that the orbital shrinkage time scale, Tghrink = ci/d^ is 
substantially longer than dynamical timescale of the secondary Tdyn,2, 



'^shrink 



1 



(8) 



where p2 is the average density of the secondary. We found that spurious oscillations are minimized for values of 
e < 0.05. As soon as the first SPH particle crosses the Lagrange point Li, we stop the relaxation process, transform 
to a space-fixed frame with the same coordinate origin and start the dynamical simulation. This moment is adopted 
as the time origin {t = 0) of the simulation, shown in Figure |2] for a 0.9 + 1.2 Mq system. In Figure [s] we show the 
potential of the initial conditions, where <l>pm is the usual Rocne potential as calculated for a point mass binary. 



^pm(r) 



GMi 
r-ri\ 



GM2 
r -r2\ 



2 (^k X r) 



(9) 



with M1/M2 and r 1/7-2 the stellar masses and centre of mass positions and the point mass Kepler frequency. The 
quantity ^ is the corresponding quantity as calculated for the actual fiuid configuration 



1 ^ 
^(r) = ^ {r) - -{uj X r) , 



(10) 



where uj is the orbital frequency of our (fiuid) binary and <I> the gravitational potential due to the fiuid. 

As we will show below in detail, the accuracy of the initial conditions has a strong impact on all major aspects of 
the simulations, prolonging the mass transfer stage by dozens of orbits in all investigated cases. 

3.2. Numerical treatment 

In this section we explore the sensitivity of the binary evolution on physical and numerical factors, including numerical 
resolution, shock heating treatment, the artificial viscosity prescription, and the initial conditions. 

3.2.1. Dependence on numerical resolution 

The duration of mass transfer and the separation at which it begins depends on the extent to which we can resolve 
the outer layers of the WD donor. In SPH, where usually the particle masses are not evolved in time, this means that 
even in a decent numerical effort with a million SPH particles the initial resolvable mass transfer rate is already a 
substantial fraction of the total system mass. The numerically resolvable mass transfer can be estimated as 



1 particle mass 



2 X 10" 



106 



-1 



orbital period ' '"^^ J \1 Mq J x 10^ cm/ s ' 

^ As the velocities vanish for synchronized systems, Coriolis forces can (and should) be ignored in this case to reach equilibrium faster. 



3/2 



ao 



\ 3/2 



Mr, 







(11) 
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Fig. 4. — Particle distribution for a test run with (upper row) and without (bottom row) shock heating. Shown are snapshots after 1, 3 
and 5 orbital periods. 

where ao is the separation between the stars at the onset of mass transfer. In order to better resolve regions of lower 
density, SPH particles of different masses could be used, but we refrain from such a possibility as it can introduce a 
substantial level of numerical noise. 

3.2.2. Comparison with the Star Crash code 
As a validation of our code we have also performed a comparison with the results obta ined with the StarCrash code. 



This code was developed originally by Rasio and Shapiro to study merging binaries rtRasio & Shapiro 1992 1994 



|T995). The major conceptual difference in their code is the use of a fast Fourier transform instead of a tree to calculate 
the gravitational forces. 

To enable a valid comparison, for this test we also use a polytropic EOS, P = Kp^^^, with F = 5/3, we evolve the 
entropic function according to 

dK^ F - 1 



dt 






■ViWij, 



(12) 



where mj is the mass of particle j, Uij the artificial viscosity tensor depending on parameters a and /3, Vij = Vi — Vj the 
velocity difference between particle i and j, and W is the smoothing kernel. We fix our artificial viscosity parameters 
in this comparison to the values recommended in their StarCrash documentation {a = 0.5, P = 1). We use the 
subroutines of the StarCrash code to set up a synchronized binary system of mass ratio q = 0.5 with 40,000 SPH 
particles, both simulations start from these initial conditions. The codes yie ld nearly indistinguisha ble results, after 
about 15 orbital periods the donor WD becomes tidally disrupted. Although Rasio & Shapiro (1995) also utilizes the 
StarCrash code, they start with slightly different initial conditions, with the two stars being closer to one another at 
t = 0. As a result, they find that the donor only survives 5 orbits, despite using the same EOS and SPH particle 
number. We suspect that they started from a slightly too small initial separation, but the very good agreement of the 
two codes for the sam e initial conditions is encouraging. In passing we note that these results are very similar to those 
D'Souza et al. (20061), which will be discussed in more detail in Section 4.4 



3.2.3. Effect of shock heating treatment 

As mentioned previously, recent grid-based simulations fMotl et al. 2002 D'Souza et al. 2006| Motl et al.|2007 ) showed 
long-lived mass transfer phases, in contrast to the earlier SPH simulations. In these calculations artincial viscosity 
was only used in the momentum equation, but not in the energy equation, which means that entropy production in 
shocks was neglected. It was suspected that at least part of the difference in the ma ss transfer duration may be due 
to the suppression or emergence of shocks/ignoring the change of entropy ( Fryer fc _D iehl 20 Q8|). 



We perform two test calculations at the example of a g = 0.5 binary, similar to D'Souza et al. (2006). We follow 
their approach and use a F = 5/3-polytrope, but in one case we evolve the entropy function according to Equation 
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Fig. 5. — Left: specific angular momentum evolution (of accretor plus disk) for the test case where we explore the effect of the (non- 
) inclusion of shock heating onto the orbital evolution. Right: evolution of the total energy dissipated, £^diss5 by the different artificial 
viscosity prescriptions. See text for details. Time is measured in units of the initial orbital period, Torbit- 



(12) and in the other we keep K constant. A mass ratio of 0.5 results in a direct impact of the accretion stream 
onto the accretor surface. In the setup that does not include shocks, the accreted matter is smoothly incorporated 
into the accretor (Figure [4j lower row), while the setup that includes shocks produces a high-entropy halo around the 
accretor (Figure [i] upper row). This high-entropy halo applies a greater torque on the binary as it has a larger effective 
lever arm than trie low-entropy halo that is produced when shocks are neglected. As a result, the matter within the 
high-entropy halo is able to remove more angular momentum from the orbit over the same number of orbital periods 
(Figure Isl) , leading to a more rapid evolution of the binary sep aration. This su pports the findings of Fryer & Diehl 
rt2008) that the large increase in orbital separation reported by [D'Souza et al.| ( |2QQ6J is in part an artifact ot shock 
suppression. 

3.2.4. Impact of the artificial viscosity prescription 

As we have shown that the entropy generation associated with shocks can affect the stability of double degenerate 
systems, we need to consider the importance of other entropy-generating mechanisms. To this end we compare 
different impl ementations of artificial viscosity. All have the "standard" form originally suggested by Monagha n fc:| 



Gingold (1983), but differ in terms of constant or non-constant parameter values and additional "limiter switches" 



• AVi: constant parameters a = 1 and (3 = 2 

• AV2: constant parameters a = 2 and (3 = 4 



AVs'. time-dependent parameters ("Morri s fc Monaghan||1997 ) , controlled as described in Rosswog et al. (2008) 
but without "Balsara switch" ([Balsara 1995|) 



• AV4: time-dependent parameters plus Balsara switch as described in Rosswog et al. (2008). 



As test case we run simulations with a 0.3 He and a 0.6 CO WD, each one modeled by 20,000 particles (T4-T7, Table 
H). Synchronized initial conditions are constructed as described in Section 3.1.2 and all tests used the Helmholtz 



OS. Artificial viscosity does indeed have an impact: in the lowest viscosity case, AV4, the system merges after 19, 
in highest- viscosity case, AV2, after 22 orbital periods. To quantify the dissipation we calculate the thermal energy 
produced by artificial viscosity before merger 



where 



diss 

dui 
~dt 



V dt J^^^ 



diss 



3 



(13) 



(14) 



is the artificial viscosity contribution to the energy equation of particle i. The quantity £^diss as a function of time 
(for t < tdisr where tdisr is the time when the donor star is tidally disrupted) is shown in Figure [s] (right panel), for 
all four artificial viscosity prescriptions. Interestingly, the prescription with the lowest average artificial dissipation 
parameters, dissipates the most energy. The reason for this is that the combination of low viscosity parameters 
and low resolution allow the SPH particles to build up a large level of velocity noise (i.e. fiuctuations around the 
local average value). This leads to an increase in the typical value of v^j, and thus more dissipation, despite the lower 
viscosity parameters and suppression of dissipation in shear fiows. 
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The increased rate of dissipation results in a more efficient spin-up of the accretor, with the angular momentum 
for this spin-up coming at the orbit's expense. This leads to a reduction in the number of orbital periods the binary 
can survive prior to merger. We perform our production runs with prescription AV4 which yields a lower limit on the 
duration of the mass transfer phase. 

3.2.5. Dependence on initial conditions 

As we have discussed, most of the earlier WD merger simulations were constructed using approximate initial condi- 
tions where the separation was estimated by comparison of the stellar radius with the Roche lobe size. To illustrate 
the importance of carefully constructed initial conditions, here we compare the evolution of a binary initialized using 
both "approximate" and "accurate" initial conditions (ICs). 




■ Accurate ICs 



■ Approximate ICs 



■ Accurate ICs 



■ Approximate ICs 




Fig. 7. — Same as Figure [6] but for the 0.6 + 0.9 Mq system. 



To construct the approximate ICs, we equate Eggleton's formula (Eggleton 1983) for the Roche radius with the 
radius R2 of the secondary as found from the SPH particle positions. I'he initial separation ao is thus given by 



Egg 



0.6(7^/^ + ln(l + 
0.49(7^/3 



(15) 



The stars, previously relaxed in isolation (i.e. without the companion's tidal field), are subsequently placed at this 
mutual separation with a spin that has the same period as the orbit of the binar y. W e refer to these ICs as "ap- 
proximate" and compare them to the carefully constructed ICs described in Section 3.1.2[ The described procedure is 
plausible, as we will see below, its approximate nature introduces a slew of artifacts. 

Two representative examples of binaries are compared, a 0.2 Mq He + 0.8 Mq CO binary and a 0.6 Mq CO + 0.9 
Mq CO binary. Each system is simulated using both the approximate (corresponding to runs Tl and T2 in Table [l] 
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Fig. 8. — Evolution of thermodynamic quantities for a 0.2 + 0.8 Mq system both with approximate (red) and accurate ICs (black). Shown 
are: total internal energy (left), maximum density (center) and temperature (right) as a function of time t in units of minutes. 
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Fig. 9. Similar to Figure[8]but for the 0.6 + 0.9 Mq system. 

and Figure [T]) and accurate (PI and P5) ICs. Some marked differences are seen between the simulations initialized 
using these two different prescriptions. The approximate prescription systematically underestimates the separation 
at which mass transfer sets in by about 14% in both cases. This is mainly because it neglects the tidal deformation 
of the donor star. As a result, binaries initialized using approximate ICs carry a deficit of angular momentum when 
compared to those evolved using the more accurate prescription. The difference in angular momentum between the 
two IC prescriptions is 11% and 13% for the 0.2 + 0.8 and 0.6 + 0.9 cases, respectively. This slight discrepancy 
severely impacts their subsequent evolution, with the duration of the mass transfer phase being underestimated by 
about two order of magnitudes. In the 0.2 + 0.8 Mq case, the approximate ICs yield two orbits (478 s) of mass transfer 
while the accurate ICs show a mass transfer for as long as 84 orbital periods (20,328 s). The difference in the 0.6 + 0.9 
Mq case is similarly dramatic (1 vs 29 orbital periods). 

The effect of varying ICs on the orbital dynamics is most clearly seen when comparing the resulting gravitational 
wave signals (Figures [6] and [7|). The approximate ICs give rise to a rapidly fading signal, while in the more accurate 
case the amplitude remains essentially constant up to the point of tidal disruption. The 0.6 + 0.9 Mq case, which 
is less stable than the 0.2 + 0.8 Mq system (Figure [T]), shows only a slightly noticeable increase in the gravitational 
wave amplitude up to disruption. In both cases, the approximate ICs overestimate the peak amplitude and frequency. 
The approximate case gives an orbital frequency of 46 mHz in the 0.6 + 0.9 Mq case and 10 mHz in the 0.2 + 0.8 Mq 
case, while the accurate ICs yield 38 and 8.4 mHz, respectively (a difference of about 20 %). This is mainly due to 
the difference in separation at final plunging (40 % difference for the 0.2 + 0.8 Mq case and 30 % for the 0.6 + 0.9 Mq 
case). 

As a result of the more abrupt artificial decrease in separation, the approximate initial conditions also give rise to 
large spurious oscillations of the central remnant core. Such oscillations are imprinted not only on the gravitational 
wave signal (see Figure [7| but are also present on the thermodynamical trajectories of the merger products. As shown 
in Figures |8] and 9) both peak densities and temperatures are overestimated when using the approximate ICs. The use 
of accurate ICs helps to better characterize the time evolution of the temperature, density, and angular momentum 
during mass transfer as well as the final morphology of the merger product, as seen in Figure p^O) Due to the larger 
initial separation for the more accurate ICs, the matter possess more total angular momentum, which results in more 
mass in the trailing arm and less mass in the disk than when starting from the approximate IC. The materi al ejected 
on eccentric or bits will fall back later onto the remnant (similar to the case of neutron star mergers; see e.g. "Rosswog] 
2007 Lee fc Ramirez- Ruiz ^^2 00 7) , the timescale of this fallback matter being seriously underestimated when using the 
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Fig. 10. — Remnant structure resulting from the merger of a 0.6 + 0.9 Mq binary system, captured 3 orbits after plunging. The left 
(right) panel shows the result for approximate (accurate) ICs. The initial angular momentum content in the system in this case differs 
by 13% between the two different IC prescriptions. This slight discrepancy severely impacts their subsequent evolution, with the duration 



of the mass transfer phase being severely underestimated by the approximate prescription, 
plunging while the accurate ICs show a mass transfer for as long as 29 orbital periods. 
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Fig. 11. — Dynamical evolution of a 0.2 He and a 0.8 Mq CO WD (at 71, 83 and 89 times the initial orbital period). The mass transfer 
in this system continues for as long as 84 orbital periods until the He WD is finally destroyed. Times are shown in the upper right corner 
of each frame in units of the initial orbital period. 

approximate ICs. 

4. WHITE DWARF MERGER SIMULATIONS WITH ACCURATE INITIAL CONDITIONS 

The use of accurate initial conditions and a correct treatment of shock heating is necessary for a reliably temporal 
evolution of the temperature, density, and angular momentum during the mass transfer phase. The evolution of two 
representative binary systems is described here in detail: the merger of a synchronized binary system with a 0.2 He 
WD and a 0.8 CO WD (PI, Table[l]) and a tidally locked binary with 0.6 + 0.9 Mq CO WDs (P5). Their evolution is 
then subsequently compared with the wide range of binary systems that we have simulated as well as with the results 
of earlier work. 

4.1. 0.2 + 0.8 Mq: A borderline unstable system near the disk stability limit 



According to the stability analysis of Marsh et al. (2004) a 0.2 He + 0.8 CO Mq binary system (marked as "PI" 
in Figure fll should still be in the "direct impact" regime where the pericentre separation of the transferred mass is 
smaller than the accretor radius, but near the stable, disk forming regime. This is clearly seen in simulation snapshots 
of the matter distribution plotted in Figure 11 In this case, we find that the rotation of the donor star remains 
synchronized with the orbital motion throughout the mass transfer phase up to the point of tidal disruption. More 
importantly, as soon as mass transfer sets in, the orbit is observed to become rapidly eccentric and the mutual stellar 
separation begins to oscillate (with about 5% amplitude changes), see Figure 12 This oscillation, as expected, also 
modulates the mass transfer rate (Figure 12). 



The average orbital separation in this system is observed to secularly increase and, as a result, the mass transfer is 
stabilized, similar to what is expected for AM CVn stars. The transfer rate stays nearly constant at a level of ~ 10^*^ 
^Edd for about 50 orbital periods and continues until the donor star is completely disrupted. This occurs after as 
many as 84 orbital or 16,930 dynamical time scales of the accretor. Before being disrupted, the donor star looses about 
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Fig. 12. — Evolution of orbital separation (in 10^ cm; left) and mass transfer rate (in units of the Eddington value; right) of the O.2+O.8M0 
binary, PI in Table [l] 
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Fig. 13. — Evolution of spin, disk and orbital angular momentum of the 0.2 + 0.8 Mq case (PI in TableQ. Jtot is the initial total angular 
momentum of the system. 

0.08 Mq. Our results are in agreement with the Marsh et al. analysis, which suggests the mass transfer phase to be 
long-lived but still unstable. 

The binary system during the mass transfer phase sheds matter through the outer Lagrange points, resulting in the 
loss of approximately 4% of the total angular momentum. As we have argued, the orbital dynamics is very sensitive 
to the loss of angular momentum which, in turn, prevents efficient mass transfer stabilization. The orbital angular 
momentum that is carried away is shared between the accreting star and the disk (Figure 13), and, due to the constant 
change in orbital separation, is observed to oscillate (Figure 13). 

The densities and temperatures of the final remnant core are shown in Figure [l4| for the 0.2 + 0.8 case. The 
highest temperatures are reached in the shearing region between the central core gina the surrounding thick disk. The 
peak temperatures of around 4 x 10^ K are reached at densities of only about 10^ g cm~^ , so that no efficient carbon 
burning is expected to occur. 

4.2. 0.6 + 0.9 Mq : An unstable binary system with direct impact mass transfer 

A 0.6 + 0.9 Mq binary system is predicted to be in the unstable, direct impact regime (P5 in Figure [T]). This is 
seen in the three dimensional density renderings plotted in Figure [Tsj which show the binary system at three different 
stages of its merger evolution. Although clearly unstable, the mass transfer continues for about 29 orbital periods 
before disruption. In contrast to the 0.2 + 0.8 Mq case, the orbit systematically shrinks from the onset of the mass 
transfer phase (Figure 16) and shows no noticeable oscillations. The depletion of the donor star in the 0.6 + 0.9 Mq and 
0.2 + 0.8 Mq binary systems is shown in Figure 17 In the 0.2 + 0.8 Mq case, the donor mass decreases very regularly 
to about 90%, followed by a phase where the mass loss speeds up dramatically and the donor becomes subsequently 
disrupted during the remaining 20 orbits. Once mass transfer begins in the 0.6 + 0.9 Mq case, the rate of mass loss 
from the donor star increases exponentially with time, leading to a swift depletion of the donor. 

The motion of the SPH particles within the Roche-type potential ^ (see Equation 10) are shown in Figure 18 As 
expected from a direct impact configuration, the angular momentum carried across is efficiently converted into the 
spin of the accreting star (Figure 19), with a small fraction being stored in the accretion torus. Similar to the 0.2 + 0.8 
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Fig. 14. — Densities and temperatures of the final remnant resulting from the merger of 0.2 + 0.8 Mq binary (PI): xy-plane (left) and 
xz-p\ane (right). Color-coded is the temperature (in units of 10^ K), the overlaid white contours refer to logp (p in g cm~^ ). The contours 
in the left panel show densities ranging from log^o p = 4 (innermost contour) to log^^Q p = 2 (outermost contour) in steps of 0.5, while the 
contours in the right panel range from log^o p = 4.1 to log^o P = 2.9 in steps of 0.3. 
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Fig. 15. — Evolution of the dynamically unstable system of two CO white dwarfs with 0.6 + 0.9 Mq. The panels show three-dimensional 
renderings of the density at 16, 27 and 30 times the initial orbital period. 
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Fig. 16. — Orbital separation (in 10^ cm; left) and mass transfer rate (in units of the Eddington value: right) for the unstable 0.6 + 0.9 
Mq system, P5 in Table [l] 

Mq binary system, the donor star is observed to remain close to synchronization during the entire duration of the 
mass transfer phase. 



The thermodynamical properties of the remnant core are displayed in Figure 20 The peak temperatures now reach 
about 5 X 10^ K at densities of about 10^ g cm~^. Under these conditions, carbon burning is observed to occur, but 
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Fig. 17. — Evolution of the donor mass normalized to its initial value for both the 0.2 + 0.8 (PI) and 0.6 + 0.9 Mq case (P5). 
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Fig. 18. — Roche contours for the system with 0.6 + 0.9 Mq components. Snapshots are taken after 3, 20 and 23 times the initial orbital 
period. 
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Fig. 19. — Dynamical evolution of the 0.6 + 0.9 Mq binary system. From left to right we have plotted the spin, disk and orbital angular 
momentum (normalized to the initial total angular momentum of the system: Jtot), respectively. 

the rate of energy injection is not large enough to be dynamicahy dominant. The structure and further evolution of 
our double degenerate merger remnant products will be discussed in more detail elsewhere. 



4.3. Summary of all other simulations 

All the systems that we have simulated in this paper belong to the regime of direct impact, unstable mass transfer. 
In all cases, we find the mass transfer to survive for many orbits. For those close to the guaranteed stable disk regime 
(see Figure |T]) the mass transfer phase extends for many tens of orbits. Systems belonging to these categories are 
represented nere by our PI (0.2 + 0.8 Mq) and P2 (0.3 + 1.1 Mq) simulations. The 0.3 + 1.1 Mq binary system shows 
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Fig. 20. — Densities and temperatures of the final remnants of the 0.6 + 0.9 Mq case (P5): xy-plane (left) and x2;-plane (right). Color- 
coded is the temperature (in units of 10^ K), the overlaid white contours refer to logp (p in g cm~^ ). The contours in the left panel show 
densities ranging from log^o p = ^ (innermost contour) to log^^Q p = 2 (outermost contour) in steps of 0.75, while the contours in the right 
panel range from log^^g P — ^-^ ^ogio P — ^-^ steps of 0.5. 

a very similar evolution to the 0.2 + 0.8 Mq system described above. The mass transfer phase is stabilized, i.e. the 
rate of mass transfer does not increases exponentially but is observed to have a nearly constant mass transfer rate for 
many orbital periods. 

The evolution of all other systems, which are far away from the limit of stable disk regime, is similar evolution to 
that of the 0.6 + 0.9 Mq binary system discussed above. In this case, a slow decrease in the separation is accompanied 
by a severe mass transfer rate increase. During the evolution, the donor star remains synchronized and the accretor is 
substantially spun up. All these systems lack stabilizing effects and a re therefore prone to K elvin- Helmholtz-tri^^ered 
He-detonations at the accretor surface that were discussed in detail in ( Guillochon et al.|2010[ ) (see Section [ 6|for further 
detai l s). Su ch systems could be possible candidates for sub-Chandrasekhar double detonation scenarios (|Fink et al. 



2007 



2010|. 

For all systems, we observe pronoun ced sp iral arms extending out of the disk which engulfs the central remnant 

this feature may be much less pronounced in 



10 



spi: 

core. As we have discussed in Section 3.2.5| and shown in Figure 
simulations that start from approximate ICs since they contain substantially less angular momentum. These cases do 
not settle into a quasi-steady-state. Instead the torques from accretor and tidal tail continue shaping the dynamical 
evolution of the disk, see e.g. Figure [TT] 

4.4. Comparison with other mass transfer calculations 

The simulations of D'Souza et al.| were the first to show that the binary system does not merge within few orbital 
time scales, but instead can survive for more than 30 orbital periods. These authors investigated the evolution of a 
polytropic system with mass ra tio q = 0.5 start ing from accurate initial conditions constructed in the framework of the 
self-consistent field technique. |D'Souza et al.| found an increase in the orbital separation and a concordant decrease 
of the mass transfer rate after about 20 orbital periods. After 32 orbits they had to stop their simulation due to the 
severe degradation of numerical conservation. At this stage, the binary system showed no sign of a pending merger. 
The authors experimented with a higher degree of contact of the same binary and consequently a higher initial mass 
transfer rate by extracting angular momentum at a rate of 1% per orbit. The merger is avoided even if the system is 
artificially driven together for three more orbits after it has reached contact. Only when driving is applied throughout 
the simulation, the donor is observed to get disrupted. This occurs after about 8 orbits. 

We found that if we reduce the angular momentum at a 0.5%/orbit rate, the disruption process is 3 times shorter, 
compared when no driving is applied. For our q = 0.5 case (0.3 + 0.6 M©, P4), for as long as we have the possibility 
to compare, we see a similar behavior to that of D'Souza et al.| in both the orbital separation and the mass transfer 
rate. At the end, the donor star is, however, disrupted after 45 orbits of mass transfer. 

The same authors also studied the evolution of a (7 = 0.4 system. They found a peak in the mass transfer rate 
after about 10 orbits followed by a rapid decrease. The authors stopped the simulation after 43 orbits, with no visible 
sign of a merger. From our calculations, P3 is the closest to their mass ratio (~ 0.41). Our P3 simulation, however, 
shows a donor disruption after 31 orbits. This difference may be due to the different shock treatment and the different 
equation of state (Helmholtz vs. 5/3 polytrope). 

5. GRAVITATIONAL WAVE SIGNALS 

The Laser Interfero meter Space Antenna (LISA) will provide th e largest observational sample of (interacting) double 
white dwarf binaries ( Farmer fc Phinney|2003 Ruiter et al.|201Q ). Such systems enter the LISA observational window 
(0.1 mHz - 100 mHz) when they reach a period ~ 5 hr. 'I'hey subsequently secularly evolve through radiation reaction 
across the LISA band and are so numerous to create a stochastic foreground. When a binary reaches a gravitational 
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Fig. 21. — Gravitational wave amplitudes of the studied double WD systems for an assumed distance of 10 kpc. The overlaid curve 
shows the LISA sensitivity for a one-year integration period that has been produced with the online sensitivity curve generator from 
http:/ /www. srl.caltech.edu/~shane/sensitivity/ 

wave frequency of a few mHz and is close enough to be resolved individually, the recorded signal shows an intrinsic 
frequency evolution over a typical one year observation time. Figure 21 shows the gravitational wave amplitudes (for 
an assumed distance of 10 kpc) for all our production runs together we the LISA sensitivity curve. Since none of the 
orbital frequencies change by more than 1 % in the last year prior to merger, we chose this as a typical integration 
time. Provided they are in the Milky Way, all of the investigated systems would be detectable by LISA. 

We show that, depending primarily on the initial conditions, the frequency evolution of the binary can be calculated 
incorrectly. Simulations that use the approximate initial conditions, as discussed in Section [3. 2. 5| severely underesti- 
mate the initial separation when mass transfer sets in. As a result, the binary only survives for a few orbits resulting 
in a rapidly fading gravitational wave signal (see e.g. Figure [t]). The accurate initial conditions, on the other hand, 
produce binary systems in which the mass transfer phase is dramatically extended by almost two orders of magnitude 
in time, resulting in a gravitational wave signal with amplitude and frequency that remain essentially constant up 
until merger. Population synthesis models have shown that these binaries are more numerous than dynamically stable 



systems for a given mass ratio (Ne lemans et al.|2QQlbl kuiter et al. 2010). As we show that these binaries can survive 
at small separation for hundreds of orbital periods (kieep in mind our mass transfer duration estimates are -due to 
finite numerical resolution- strict lower limits on the durations realized in nature), their associated gravitational wave 
signal should be included when calculating the LISA gravitational wave foreground (although expected to below LISA's 
sensitivity at these high frequencies). 



Hydro Point masses Hydro Point masses 




Fig. 22. — 0.2 He and 0.8 Mq CO WDs: comparison of the gravitational wave amplitude /i+ of the simulation with the point mass limit, 
r is the distance to the observer. 



With the use of accurate initial conditions as well as realistic shock heating treatment, the frequency evolution of the 
binary can thus be explored in detail. Depending primarily on the mass ratio, as we have argued here, the resulting 
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Fig. 23. — CO WDs with 0.6 + 0.9 Mq\ comparison of the gravitational wave amphtude /i+ of the fuh hydrodynamic simulation with 
the point mass limit, r is the distance to the observer and t is the time measured in units of the initial orbital period. 

gravitational wave signal can be significantly different from that predicted in the point mass limit, particularly for 
marginally unstable systems. In such systems, the gravitational wave signal carries the imprint of the orbital oscillation, 
as illustrated by the evolution of the marginally unstable 0.2 + 0.8 system in Figure 22 



For systems that are far from the stable disk limit, we find the point mass approximation to provide a better 
description of our simulated gravitational wave signals. This is mainly because the oscillations in the orbital separation 
are much less pronounced for these systems. Contrary to the 0.2 + 0.8 case, the gravitational wave amplitude 
derived for the 0.6+0.9 Mq binary system remains close to the point mass prediction up to almost complete disruption, 
as seen in Figure |23j 

6. SURFACE DETONATION CRITERIA IN WHITE DWARF BINARIES 

During the final phase of mass transfer prior to merger, the rate of mass exchanged approaches the mass of the 
donor divided by the orbital timescale. 




O.IMqS 



(16) 



where \i\ is the reduced mass of the accretor. As shown in Guillochon et al. (2010), Kelvin-Helmholtz instabilities 
develop within the accretion stream once the accretor has built up a torus of material acquired from the donor, and 
this leads to the formation of dense knots of material that can periodically strike the accretor's surface. This striking 
action can lead to the initiation of a detonation within the torus about the accretor, assuming that matter acquired 
from the donor mostly consists of helium. 

The question of whether knots are capable of detonating the accretor's helium layer depends on the ram pressure 
applied by individual knots, and by the initial conditions of the helium torus at the time of compression. From the SPH 
simulations presented in this paper we have an accurate record of the time-evolution of the accretion stream's density 
at LI and the conditions at the base of the helium torus. As shown in our FLASH simulations (Figure 25), a standing 



shock develops between the stream and the torus, equalizing the pressure at the stream-torus interface, and reducing 
the velocity of the fluid flow in the torus relative to the accretion stream to subsonic velocities. As the interaction 
is subsonic and no heat transport processes can operate effectively on a dynamical timescale, the accretion stream 
smoothly evolves relative to the background pressure, and thus maintains constant entropy. The exterior pressure 
experienced by the accretion stream is given by the pressure within the torus. 

Given the initial density pLi and temperature Tli of the stream material as it leaves LI, the pressure at the base 
of the helium torus Ptor, and an equation of state that parameterizes the entropy s and pressure P as functions of 
density and temperature, we can determine the density of the stream when it reaches the base of the helium torus by 
solving the following system of equations for pknot and T^not, 



5(PL1,Pli) = 5(pknot,Pknot) 
P(p 

tor? Ptor 



(17) 
(18) 



where "tor" refers to the base of the helium torus and "knot" refers to the stream post-compression. For all of the 
SPH simulations presented in this paper, the maximum increase in density achieved within the stream as it falls from 
LI to the accretor's surface is > 10. The small initial entropy of the stream (Tli ^ 10^ K) and minimal shock-heating 
means that the final temperature reached within the stream is substantially smaller than the temperature required for 
dynamical burning, and thus burning is not expected to be important within the stream itself. 

With our estimate for P knot: we ca n estimate the increase in density and temperature of the torus after it is compressed 
by the stream. As in , Guillochon et al.^(2010J , we solve for the trajectory of a test particle released at LI with initial 
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Fig. 24. — Side-by-side comparison of the fluid internal energy pe for one of the SPH calculations featured in this work (Gl) and one of the 
FLASH calculations of ( Guillochon et al.|2010, run Ba) for a 0.45 + 0.9 Mq binary. The SPH calculation fully resolves both the donor and 
the accretor, which allows tor the question ot binary stability to be addressed self-consistently. However, the SPH calculation has limited 
resolution and is unable to sufliciently resolve the stream-torus interface, artiflcially suppressing the instabilities and the standing shocks 
that are observed to develop in the FLASH calculation. As a result of not resolving the standing shocks properly, the SPH calculation 
shows somewhat less heating in the torus than the FLASH calculation. On the other hand, the stream boundary condition used in the 
FLASH simulation is idealized to lie on the Roche surface of the donor, with a simplifled initial velocity proflle normalized to the sound 
speed of the donor, which does not perfectly represent the actual flow of material leaving LI. 




Fig. 25. — Ratio of the triple-alpha burning timescale tsq, to the dynamical timescale Tdyn for different combinations of donor and accretor 
masses. The flgure shows the same stability criteria as Figure ^ sans labels for clarity. The colored region shows the ratio of the two 
timescales, with the white contours showing where tsq, = O.lrdyn, tsq, = Tdyn, and tsq, = lOrdyn- The dotted line shows the convex hull of 
all simulations performed in which the donor can possess a signiflcant amount of helium in its outer layers. For the two FLASH simulations 
performed in Guillochon et al. ( 2010) in which a surface detonation was observed (Gl and G3), the ratio of timescales is of order unity, 
whereas the simulation m which no detonation was observed (G2), tsq, is signiflcantly larger than Tdyn- 

velocity given by the sound speed of the donor to estimate v^^ the component of the velocity perpendicular to the 
accretor's surface. Given the density of the stream pknot, we can determine the ram pressure of the stream 



Piam — Pkn 



(19) 



We can then use Pram 
state of the torus 



and solve a set equations similar to Equations (17) and (18) to determine the post-compression 



<^(Ptor 7 Ptor) — "^(Pcomp? Pcomp) 
-Pram ~ P(Pcomp7 Pc( 



(20) 

"^comp? J- comp (21) 

where "comp" refers to the conditions in the helium torus after compression. This enables us to calculate the triple- 
alpha burning timescale t^cx and compare it to the dynamical timescale Tdyn at the accretor's surface to determine if 
runaway burning can be achieved, leading in turn to a detonation. 

In Figure 25 we show the curves for which Tdyn = tsq, and lOrdyn = T^^a within the co nvex hull of the SPH sim ulations 
presented inmis paper. For the two simulations which exhibited surface detonations in Guillochon et al. (2010) the two 
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timescales are comparable, whereas tdyn is far shorter than r^a for the simulation that showed no surface detonation. 
The SPH simulations show that the conditions for dynamical burning being triggered on the surface of the accretor 
are achieved for accretor masses larger than ~ O.SM©, with a cutoff at low donor masses as the circularization radius 
of the accretion stream begins to exceed the radius of the a ccretor. Such sys t ems c ould be possible candidates for 



sub-Chandrasekhar double detonation scenarios, as argued in Guillochon et al. (2010) 



7. SUMMARY 

The onset of (numerically resolvable) mass transfer in double degenerate binary systems and its impact on the orbital 
evolution forms the main theme of our paper. We find for all investigated cases that mass transfer ensues for at least 
dozens of orbits, even for systems which according to the analysis of Mar shetahj ( 2004 ) are guaranteed to be unstable. 
Due to the finite numerical resolution, the duration of the mass transfer phase calculated here should be considered 
as a strict lower limit. 

Part of our motivation has been to settle a previous debate concerning the duration of mass transfer in such systems. 
Almost all of the earlier SPH simulations had found quick disruption of the donor star within a few orbital periods 
after the onset of mass transfer, while recent grid-based simulations ( jP'Souza et al.||2006 ) found a longer lived mass 
transfer phase. Using the same assumptions and initial conditions we have reproduced essentially all previous results. 
In addition we have compared our results to those obtained by using the StarCrash code developed by Faber and Rasio 
and found excellent agreement when using the same initial conditions and equation of state. 

Concerning the mass transfer duration debate, we find that the reported discrepancies can be accounted for by the 
following two effects. The first is the self-consistent inclusion of entropy production in shocks, which we find to have 
a major impact on the orbital dynamics. In direct-impact systems, we find the realistic inclusion of shock heating to 
produce an extended, high-entropy halo around the accretor. When ignored, a smooth incorporation of the transferred 
material into the stellar surface is observed. Since substantial amount of angular momentum can be stored in the 
shocked material (which is provided at the expense of the orbital angular momentum), the rea hstic inclusion of shoc ks 



yields shorter lived mass transfer phases. These results confirm the earlier concerns raised by Fryer & Diehl (2008). 

The second and even more important reason is attributed to the initialization procedure ot the binary m the 
simulations. The results of simulations constructed using approximate corotating initial conditions have been compared 
with those using carefully constructed binaries and are found to be markedly different. Simulations that use the 
approximate initial conditions underestimate the initial separation at which mass transfer sets in by about 15%, 
which yields a binary that only survives for a few orbits and thus a rapidly fading gravitational wave signal. This 
is because once mass transfer sets in, the orbital evolution is no longer driven by angular momentum losses due 
to gravitational wave emission but instead is dominated by the redistribution of angular momentum in the system. 
The use of approximate initial conditions systematically underestimates the number of orbits and overestimates the 
gravitational wave peak frequencies and amplitudes by about 20 %. In contrast, the accurate initial conditions produce 
a binary system in which the mass transfer phase is extended by almost two orders of magnitude in time, resulting 
in a gravitational wave signal with amplitude and frequency that remain essentially constant up until merger. The 
mass transfer also shows a unique oscillatory signal that is most pronounced for binary systems that are near the disk 
forming limit. 

The use of accurate initial conditions and a correct treatment of shock heating allows us to better characterize the time 
evolution of the temperature, density, and angular momentum, which are important when considering thermonuclear 
events that may occur during the mass transfer phase and/or after merger. Approximate initial conditions, for example, 
tend to significantly overestimate the densities and temperatures of the final merger remnant. This fact may have 
adverse consequences for possible type la candidate systems. Our treatment allows us to also accurately identify when 
surface detonations may occur in the lead-up to the merger. The results of this study have been used as boundary 
conditions for FLASH simulations that focus o n the fate of the low-density accreted material and that are described 
in a companion paper ("Guillochon et al. 2010). This combined study has shown that contrary to earlier beliefs, the 



mass transfer can be very eventful. We have shown that for the case of unstable, direct impact He accretion onto the 
accretor a hot (> 5 x 10^ K) helium torus builds up in which the incoming accretion stream triggers Kelvin-Helmholtz 
instabilities. Such Kelvin-Helmholtz "knots" repeatedly impact the accretor's helium surface and, once the triple-alpha 
time scale is shorter than the local dynamical time scale, can trigger violent surface explosions. This explosion launches 
shocks into the accretor's interior that can subsequently ignite the carbon oxygen core. This could possibly be the long 
sought-after explosion mechanism how to trigger supernova explosion in a double degenerate scenario which, together 
with the promising rates, could make double degenerate mergers a viable (sub-luminous) type la progenitor system. 
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